namespace Likelihoods{

class Poisson : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Poisson(double *Y, mat X, Distribution::Distribution *S,mat s) : Density::GeomBridge(Y, X, S){ 
			_s=s;
			_is=inv(s);
			//d= new boost::math::normal_distribution<>(0,1); 
  }
		~Poisson()
		{
			//delete d;
		}
		inline double Likelihood(mat theta1)
		{
			mat theta=theta1.t();
			double L=0;
			double sum=0;
			int m=theta.n_cols;
			//Xb
			int n=Get_n();
			mat X=Get_X();
			double *Y=Get_Y();
			for(int i=0;i<n;i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= X.row(i);
				sum=dot(foo,theta);
				L+=-exp(sum)+Y[i]*sum-logfact(Y[i]);
			}
			//cout << L+12200;
			return L+12000;

		}
		double Prior(mat theta){
			double res=0;
			mat foo=theta.t();
			int d=theta.n_elem;
			mat bar=foo*_is*foo.t();
			//cout << "//" << det(_s);
			double de=det(_s);
			if(de>100000){de=100000;}
			res=-d*0.5*(log(2*PI))-0.5*log(de)-0.5*as_scalar(bar);
			return res;
		}
		mat GradLik(mat theta)
		{
			mat X=Get_X();
			double *Y=Get_Y();
			int p=X.n_cols;
			mat sum(p,1);
			sum.fill(0);
			int n=X.n_rows;
			double xb=0;
			for(int i=0;i<n;i++)
			{
				xb=dot(X(i,span::all),theta);
				
				double gpxb=Phi(xb);
                                double pxb=exp(-xb);
				
				//cout << "foo1 "<<foo1;
                                //cout << "foo2 "<<foo2;
				sum+=-(exp(xb)-Y[i])*X(i,span::all);
			}
			sum+=-as_scalar(arma::sum(theta,0))/_s(0,0);
					
			return sum;

		}
	        mat Get_s(void){return _s;}	
	private:
		mat _s;
		mat _is;
		//boost::math::normal_distribution<> *d;
	/*private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;*/
};

} // end namespace

